Sequestration of noble gases in giant planet interiors 
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The Galileo probe showed that Jupiter's atmosphere is severely depleted in neon compared to protosolar 
values. We show, via ab initio simulations of the partitioning of neon between hydrogen and helium 
phases, that the observed depletion can be explained by the sequestration of neon into helium-rich 
droplets within the postulated hydrogen-helium immiscibility layer of the planets interior. We also 
demonstrate that this mechanism will not affect argon, explaining the observed lack of depletion 
of this gas. This provides strong indirect evidence for hydrogen-helium immiscibility in Jupiter. 



Jupiter is the most extensively probed and best under- 
stood of the giant planets, but many questions regarding 
its detailed composition, formation, and interior struc- 
ture remain unanswered. One issue of major importance 
to structural models is the question of whether hydrogen 
and helium mix homogeneously throughout the planet 
or whether a layer of hydrogen-helium immiscibility ex- 
ists deep within the interior In the immiscibility 
layer, helium would form dense droplets which would rain 
down into the deeper interior and redissolve, resulting in 
a gradual and ongoing transfer of helium from regions 
above the immiscibility layer to regions below. Such a 
layer almost certainly exists in Saturn, as evident from 
the observed depletion of helium from its upper atmo- 
sphere (compared to protosolar values) and the appar- 
ent excess luminosity of the planet [fjj. For the hotter 
interior of Jupiter the case is less clear since there is 
no measurable excess luminosity and the observed he- 
lium depletion from the upper atmosphere is quite small 
(0.234 by mass compared to 0.274 in the protosolar neb- 
ula [5rtZ|). Theoretical attempts to determine the pres- 
sure/temperature range in which H and He are immisci- 
ble u sing successively more sophisticated levels of theory 
[H, 0-0] have produced quite different results, however 
recent work [l3l H3 | provides a hydrogen-helium immis- 
cibility line which is very close to the Jupiter isentrope 
in the 100-300 GPa region. 

The strongest evidence for H-He immiscibility may in 
fact come from the depletion of neon. Jupiter's upper at- 
mosphere was found by the Galileo entry probe to be ex- 
tremely deficient in neon: although neon makes up 1/ 600 
by mass in the solar system it comprises only 1/6000 
by mass in Jupiter's upper atmosphere Prior to 

these measurements, Roulston and Stevenson [lj| pro- 
posed that hydrogen-helium immiscibility could lead to 
neon depletion on the assumption that neon would pref- 
erentially dissolve in the helium-rich phase in the immis- 
cibility layer. This would lead to Jupiter's neon content 
being gradually carried down within the helium droplets 
and concentrating in the deep interior. When hot fluids 
enriched in He and Ne are converted upwards they are 
subjected to differentiation again. There is, however, a 
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FIG. 1: (left) Schematic depiction of the interior of a gas gi- 
ant (e.g. Jupiter or Saturn) with a layer of H-He immiscibil- 
ity. Helium-rich droplets form within the immiscibility layer 
and rain downwards, leading to a slow increase in the helium 
concentration in the deep interior. Neon is absorbed within 
the droplets and carried out of the upper atmosphere. (right) 
P/T curves for Jupiter and Saturn combined with the loca- 
tion of the H-He immiscibility region determined in the work 
of Morales et al 114] assuming an overall He atomic molar 
concentration of 0.0847. 



lack of direct experimental evidence for whether Ne will 
indeed preferentially dissolve in the helium phase as pro- 
posed. The pressure-temperature conditions correspond- 
ing to phase separation can currently only be attained in 
shock-wave experiments lasting only tens of nanoseconds 
16] and which are hence ill-suited to study questions of 
phase separation and partitioning. Is it also not known 
why the chemically similar noble gas argon is not seen 
to be depleted but instead is present at slightly above- 
solar concentrations (« 1.6 times solar Q) comparable 
to most other detectable trace heavy elements in Jupiter, 
and whether this indicates that the depletion mechanism 
acting upon neon does not act upon argon or whether it 
indicates a very high initial argon concentration. In or- 
der to resolve these issues, here we present ab initio free 
energy calculations with a view to determining the solu- 
bility behaviour of neon and argon in H and He phases at 
pressures corresponding to the postulated immiscibility 
region. 

The distribution of a trace species such as Ne or Ar 
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between coexisting phases is dependent upon the Gibbs 
free energy of transfer, AGtt, being the change in G when 
an atom of the trace species is moved from one phase to 
the other at constant P and T, in this case 



AGty = G(He + Ne) + G(H)-G(H + Ne)-G(He). (1) 

Here we compute AG-it for neon and argon in pure H 
and He within the density functional theory molecular 
dynamics (DFT-MD) framework, within a temperature 
and pressure range of 100-300 GPa and 3000 to 7000 K. 
Determination of free energies from MD is a nontrivial 
problem for which a number of methods have been devel- 
oped. We use a two-step coupling constant integration 
(CCI) approach |17l| similar to that recently applied by 
Morales et al. [14j , in which the Gibbs free energy of the 
DFT system is determined by adiabatically transforming 
the system in two steps: (a) from the DFT system to a 
system of atoms governed by a classical potential and (b) 
from the classical system to a noninteracting gas. 

The CCI method provides a general scheme for com- 
puting the Helmholtz free energy difference between two 
systems governed by potential energy functions U\ and 
U2. Constructing artificial system with potential U\ = 
(1 — X)Ui + XU2, the free energy of system 2 may be 
expressed as: 



F x + / d\(U 2 - 0i} A , 



(2) 



where at each integration point, the average is taken 
over a sample of configurations obtained in the U\ sys- 
tem. Since the difference in Gibbs free energy between 
a system at two different pressures can be found by the 
thermodynamic integration G(P2) — G(Pi) = Jp^ VdP, 
we performed all CCI calculations at pressures close to 
200 GPa and then integrated the equation of state for 
each system outwards to obtain G values at pressures 
from 100 to 300 GPa. 

The first part of the CCI was the integration from the 
noninteracting system to the classical system. The clas- 
sical potential used was a pair potential of a modified 
Yukawa form 1141: 



U(r) 



-br 



e -HL-r) 

~(L~r) 



J>L/2 



(3) 



for r < L/2 and zero otherwise. We set L — 9.749 
a.u., then fitted a and b to the g(r) functions of the DFT 
systems at high pressure 2J|. The integration from the 
noninteracting to classical system used sixteen A values, 
spaced more closely in the small- A region where the com- 
puted classical energy varies rapidly. We checked po- 
tentials obtained by a force-matching method and found 



that the numerical results achieved for the free energy of 
the final system were not altered by the different poten- 
tial. 

The second part was the integration from the classi- 
cal forces to the DFT forces. This was the most time- 
consuming step and required a series of DFT-MD runs. 
We found that the variation in (C/dft — ^classical) was suf- 
ficiently close to linear in A to allow a fit with only three A 
points to be used - checks against calculations with five 
A points resulted in discrepancies smaller than 0.1 eV. 
All DFT simulations were performed using the VASP 
code [18|. We used 128 H atoms or 64 He atoms with 



or without a single Ne or Ar atom inserted per cell. We 
used pseudopotentials of the projector-augmented wave- 
function type [lj} , the exchange-correlation functional of 
Perdew, Burke and Ernzerhof [20J, a cutoff of 1000 eV 
and eight fc-points in the Monkhorst-Pack grid. All MD 
simulations used a timestep of 0.4 fs and were run for 
5000 timesteps with the first 1000 steps discarded for 
equilibration. 

The total Gibbs free energy computed for each sys- 
tem (H, H+Ne,H+Ar,He,He+Ne,He+Ar) is a sum of five 
terms: 



G(Pi) = -Fideal(Vb) + APidcal->classical(Vb) + 

AF c i assica i_>DFT(^o) + J fWo+ / VdP, (4) 

JPo 

where Po,V$ are the pressure and volume at which the 
F values were computed, and P\ is the pressure of in- 
terest. The CCI procedure was undertaken at pressures 
within 1% of 200 GPa and VdP integration was used to 
correct the values back to 200 GPa exactly. Pressure- 
volume curves were obtained from a series of five MD 
simulations on each system at pressures spaced from 100 
to 300 GPa, and by fitting the resulting data points with 
a piecewise power law fit. 

In order to validate this method, we also performed 
simulations via an alternative free energy calculation 
method based on the particle insertion formalism of 
Widom (2l| . Using only the V point for Brillouin Zone 
sampling, we computed the free energies associated with 
the insertion of Ne and Ar into He and H cells at volumes 
corresponding to a Wigner-Seitz radius for the electrons 
of 2.4 bohr radii, then integrated along isotherms to ob- 
tain values of AGt r which were then compared with CCI- 
computed values. The results were found to agree within 
the relevant error bars. Since the CCI method is more 
computationally efficient we applied it for the computa- 
tion of the final, eight fc-point results. We also estimated 
the quantum correction to the classical free energy result- 
ing from the fluctuations around the classical trajectories 
of the nuclei. Using the first term of a Wigner-Kirkwood 
expansion in 7t(25|. we estimate the free energy correction 
at 5000 K and 200 GPa to be of the order of 0.01 eV per 
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T P 
(K) (GPa) 


AGtt (eV), CCI 
Ne Ar 


AG Tr (eV), IM 
Ne Ar 


200 3000 
200 5000 
200 7000 


-2.45(33) 4.88(33) 
-2.36(46) 5.08(45) 
-2.42(63) 5.59(66) 


-2.01(37) -0.10(31) 
-1.17(41) 0.41(36) 
-1.14(48) -0.38(43) 


100 5000 
300 5000 


-1.61(47) 4.73(47) 
-2.78(46) 5.65(46) 





TABLE I: Computed CCI AGir values for neon and argon 
in pure hydrogen vs pure helium phases as a function of tem- 
perature and pressure. A negative sign indicates preference 
for solubility in the helium phase. Also shown are the AGTr 
values computed using the ideal mixing (IM) approximation, 
that is, with G = U + PV — TSidcai where Sidcai is the entropy 
of the ideal gas. The resulting Gti values are significantly 
smaller than those computed with the CCI method. 



atom or less, and can consequently be neglected. 
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FIG. 2: At left, the difference in volume A Virus between the 
pure hydrogen/helium cells and the cells with a single Ne/Ar 
atom added isobarically at 5000 K. At right, the computed 
difference in free energy AGi ns between the pure H/He cells 
and the cells with a single Ne/Ar atom added for pressures 
between 100 and 300 GPa at 5000 K. 



The computed values of AGT r for neon and argon are 
given in Table fl] For neon at 200 GPa we find AGxr val- 
ues of approximately —2.4 eV with temperature variation 
from 3000 K to 7000 K producing only a small variation 
in AGTr- The negative sign here indicates a preference 
for solubility in helium. Thermodynamic integration of 
the pressure- volume curves as shown in Figure 2(a) from 
200 GPa shows that the helium preference becomes more 
pronounced with increasing pressure. For argon, we see 
contrasting behaviour: at 200 GPa and 5000K we have 
a AGTr of +5.1 ± 0.45 eV with the positive sign now in- 
dicating a preference for the hydrogen phase. The mag- 
nitude of the preference for hydrogen solubility increases 
somewhat with both temperature and pressure. 

Following the work of Roulston and Stevenson [l5| , we 



expect the rate at which neon is removed from the upper 
envelope to be related to the loss rate of helium by 



dXj^e ( AGtiA dXft e 

= Anc exp f 



dt 



V k B T 



dt 



(5) 



This implies that the observed depletion of neon will 
be approximately given by 



log 



(6) 



where Xq.q and Xq^\ refer to the original (protosolar) 
and present-day molar atomic concentrations of species 
Q in the upper Jovian atmosphere, respectively. Based 
on the measurements of Von Zahn et al. [6] for the cur- 
rent helium concentration and the estimate of Lodders 
for the protosolar concentration, we find a helium de- 
pletion Xjj e — X^ c value of approximately 1.2%. Com- 
bining this with values of AGti of —2.35 ± 0.45 eV for 
neon partitioning, we obtain the relationship between T 
and neon depletion shown in Fig. [3] The observed value 
of approximately 0.1 for the neon depletion ratio corre- 
sponds to T values of between approximately 4000 K and 
6000 K. This is consistent with the expected temperature 
of the immsicibility region. The computed value of AGti- 
is thus consistent with the assumption that the observed 
depletions of both helium and neon are due entirely to he- 
lium rain within the hydrogen-helium immiscibility layer. 
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FIG. 3: Relationship between the temperature of the immis- 
cibility region and the change in helium concentration which 
would be required to produce the observed 90% depletion of 
neon, for AGty values of -2.35 ± 0.45 eV. The observed value 
of 1.2% for A^He is marked with a line. 

For argon, the positive value of AGti- implies that Ar 
will be almost entirely excluded from the He phase. Since 
the helium phase remains only a very small portion of 
the planet this will lead only to a miniscule enhancement 
in the argon content of the upper atmosphere. This im- 
plies that the measured concentration, approximately 1.6 
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FIG. 4: Pair correlation functions g(r — i? so i) for distances 
between solvent (H,He) an d solute (Ne,Ar) atoms at 200 GPa 
and 5000 K. The curves are shifted by the effective radius i? so i 
of the solvent atom in each case, determined from the point 
where the solvent-solvent g(r) crosses 0.5. i? so i is 0.37 A for 
hydrogen and 0.58 A for helium. 



times the solar value [B[ , should be close to the true argon 
concentration of the planet as a whole. 

The difference in solubility behaviour between neon 
and argon invites further examination. The difference 
in the free energies of insertion is governed primarily by 
the volume change AVins associated with the insertion at 
constant pressure of the noble gas atom into the pure- 
solvent cell. As shown in Fig. 2(a) the effective volume 
of neon is larger by 0.86 A 3 at 200 GPa and 5000 K in 
hydrogen than in helium, while argon shows the opposite 

! 3 

trend, being larger by 0.73 A in helium than hydrogen. 
In Fig. 2] we plot the pair correlation function g(r) for 
the solvent atoms surrounding each species of noble gas 
atom. The g(r) curves have been shifted by i? S oi, the ef- 
fective radius of the solute atom derived from the point 
at which the g(r) for H-H or He-He crosses 0.5. There 
is a clear difference in the exclusion behaviour of neon 
and argon, with the small-distance tail of the H-Ar curve 
allowing a closer effective approach than in helium, in 
contrast to of the H-Ne and He-Ne curves where helium 
approaches more closely. 

As a possible interpretation, wc note that in the 
P,T range of interest the H atom is essentially ionized 
22J, |23j while the He atom retains its electrons. A he- 
lium atom thus is repelled from the Ne/Ar atom by the 
electron-electron interactions dominated by Pauli exclu- 
sion, whereas hydrogen atoms may more easily penetrate 
the outer shells and are repelled primarily by core-core 
repulsion. The Ar atom's additional electron shell thus 
gives it a larger effective volume to exclude the helium 
atom, but much less so for the hydrogen. If this model 



is correct then we would expect the noble gases krypton 
and xenon to likewise exhibit a preference for hydrogen 
solubility, a behaviour consistent with their apparent ob- 
served lack of depletion in the upper atmosphere [j| . 

In this work we have considered only pure helium and 
pure hydrogen phases. In practice the helium phase will 
have very little hydrogen, but the hydrogen-dominant 
phase still contains some helium [l4[ , however we do not 
expect this to qualitatively change the results. Another 
limitation not considered in this study is whether the par- 
titioning coefficient changes as the neon concentration in 
helium increases; it should be noted that the required 
molar concentration of neon in the pure-helium phase 
will be quite large. We also cannot exclude, based on 
this study, the possibility that neon forms its own phase, 
however we consider this unlikely due to the small initial 
Ne concentration. 

These results strongly support the existence of 
hydrogen-helium phase separation in Jupiter as an ex- 
planation for the observed Ne depletion. We have also 
shown that argon will be preferentially excluded from a 
helium droplets explaining the observed lack of depletion 
of this element. Further work to more accurately deter- 
mine the location of the hydrogen- helium immiscibility 
line at low pressure and low temperature would allow us 
to make a quantitative estimate of the neon concentration 
in Saturn to be tested by future missions. Furthermore, 
neon may be added as a tracer in laboratory experiments 
to detect the phase separation of H-He mixtures because 
neon scatters X-rays more strongly. 
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